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Abstract. The most interesting step of condensation is the cluster formation up to the critical size. 
In a closed system, this is an instationary process, as the vapour is depleted by the emerging liquid 
phase. This imposes a limitation on direct molecular dynamics (MD) simulation of nucleation 
by affecting the properties of the vapour to a significant extent so that the nucleation rate varies 
over simulation time. Grand canonical MD with McDonald's daemon is discussed in the present 
contribution and applied for sampling both nucleation kinetics and steady-state properties of a 
supersaturated vapour. 

The idea behind that approach is to simulate the production of clusters up to a given size for a 
specified supersaturation. In that way, nucleation is studied by a steady-state simulation. A series 
of simulations is conducted for the truncated and shifted Lennard-Jones fluid which accurately 
describes the fluid phase coexistence of noble gases and methane. The classical nucleation theory is 
found to overestimate the free energy of cluster formation and to deviate by two orders of magnitude 
from the nucleation rate below the triple point at high supersaturations. 
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INTRODUCTION 

The key properties of nucleation processes in a supersaturated vapours are the height 
A£2* of the free energy barrier that must be overcome to form stable clusters and the 
nucleation rate J? that indicates how many macroscopic droplets emerge in a given 
volume per time. The most widespread approach for calculating these quantities is 
the classical nucleation theory (CNT) [1], which has significant shortcomings, e.g., it 
overestimates the free energy of cluster formation [2, 3]. An important problem of CNT 
in case of vapour to liquid nucleation is that the underlying basic assumptions for the 
liquid do not apply to nanoscopic clusters [4-6]. 

Molecular simulation permits the investigation of nanoscopic surface effects and the 
stability of supersaturated states from first principles, using effective pair potentials. For 
instance, the spinodal line can be detected with Monte Carlo (MC) [7] simulation meth- 
ods; in experiments, it can only be approximated as it is impossible to discriminate an 
unstable state from a metastable state where AQ* is low. Equilibria [8] and vapourization 
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processes [9, 10] of single clusters can also be simulated to obtain the surface tension 
as well as heat and mass transfer properties of strongly curved interfaces. Moreover, 
molecular dynamics (MD) [11-13] and MC [14] simulation of supersaturated systems 
with a large number of particles are useful for the study of very fast nucleation pro- 
cesses, whereas lower nucleation rates can be calculated by transition path sampling 
based methods [15, 16]. 

Equilibrium simulations fail to reproduce kinetic properties of nucleation processes 
such as the overheating of growing clusters due to latent heat. On the other hand, direct 
MD simulation of nucleation, where cluster formation is observed directly in a near- 
spinodal supersaturated vapour, has its limits: if nucleation occurs too fast, it affects 
the properties of the vapour to a significant extent so that the nucleation rate obtained 
according to the method of Yasuoka and Matsumoto [11] and other properties of the 
system vary over simulation time [17]. In the present work, nucleation is studied as 
a steady-state process by combining grand canonical MD (GCMD) and McDonald's 
daemon [18, 19], an 'intelligent being' that eliminates large droplets from the system. 



SIMULATION METHOD 

Supersaturated states can be characterized in terms of the difference between the chem- 
ical potential /i of the vapour and the saturated chemical potential fi a (T). The chemical 
potential of the vapour can be regulated by simulating the grand canonical ensemble with 
GCMD: alternating with canonical ensemble MD steps, particles are inserted into and 
deleted from the system probabilistically, with the usual grand canonical acceptance cri- 
terion [20]. For a test insertion, random coordinates are chosen for an additional particle, 
and for a test deletion, a random particle is removed from the system. The potential en- 
ergy difference AY due to the test action is determined and compared with the chemical 
potential. The acceptance probability for insertions is 



& = min ( 1 , exp 
while for deletions it is 

& = min I 1 , exp 
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wherein A is the thermal wavelength. Of course, care must be taken that the momentum 
of the inserted particles is consistent with the simulated ensemble and does not introduce 
any artifical velocity gradients. The MD integration time step was At = 0.00404 in 
reduced time units, i.e., o(m/ e) 1//2 , wherein £ is the energy parameter of the fluid model 
and m is the mass of a particle. The number of test actions per simulation time step 
was chosen between 10~ 6 and 10~ 3 N, a value which was occasionally decreased after 
equilibration if very low nucleation rates were observed. 

Molecular simulation of nucleation has to rely on a cluster criterion to distinguish the 
emerging liquid from the surrounding supersaturated vapour [21]. In the present case, the 
Stillinger criterion [22] was used to define the liquid phase and clusters were determined 



as biconnected components. Whenever a cluster exceeded the specified threshold size 0, 
an intervention of McDonald's daemon removed it from the system, leaving a vacuum 
behind [18, 19]. 



NUCLEATION THEORY 

The free energy of cluster formation is the same for the grand canonical and the 
isothermal-isobaric ensemble [23]. At specified values of the chemical potential fi of 
the supersaturated vapour, the total system volume V and the temperature T, it is related 
to the surface energy r\ by [24] 

AQ V = ( p - pe )dVi+ U^l (ft- H)d V, (3) 
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where v is the number of particles in the cluster, p is the supersaturated vapour pressure, 
Vf (v) is the volume and ^"(v) the surface area of a cluster containing v particles. Note 
that ji( as well as pi are the chemical potential and the pressure of the liquid phase at the 
conditions prevailing inside the cluster. In CNT, it is assumed that the bulk liquid density 
at saturation p' and the density of a nanoscopic cluster are the same and all clusters are 

treated as spheres, i.e., pi = p' and &{y) = ^.(v) = (6^KV / p') 2 ^ . Accordingly, the 
chemical potential of the liquid inside the nucleus is approximated by 

w ^(r) + r^^( T)+ ^, (4) 
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and the cluster surface tension f= (dr\/d^) by the surface tension y of the planar 
vapour-liquid interface, leading to [25, 26] 



dn = 




dv. (5) 



The free energy of formation has a maximum AQ * which lies at the size v* of the critical 
nucleus. Including the Zel'dovic factor fz and the thermal non- accomodation factor /at 
of Feder et al. [1], the nucleation rate is 

J =/ 4 r/zyexp(-/3A^)^V(v*), (6) 

where ~N\ is the number of vapour molecules in the system and h is the Planck constant. 

Instead of using the surface tension of the planar interface, Laaksonen, Ford, and 
Kulmala (LFK) [27] proposed an expression equivalent to 



J yd^ = y&(v) ^l + aiv~ 1/3 + a 2 v~ 2/3 J 



(7) 



The two parameters ot\ and «2 are determined from the assumption that almost all 
particles are arranged either as monomers or as dimers and that the Fisher [28] equation 



of state correctly relates p/T to the number of monomers and clusters present per 
volume. Effectively, LFK theory modifies CNT only by the introduction of the parameter 
Ct\, since OC2 cancels out for all free energy differences if the usual assumption & ~ v 2 / 3 
is applied. 

The Hale scaling law (HSL) is based on a different approach [29]. In agreement with 
experimental data on nucleation of water and toluene [29], it predicts 
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with a proportionality constant depending only on properties of the critical point. 

In the present work, these theories are evaluated using Gibbs-Duhem integration over 
the metastable part of the vapour pressure isotherm collected by canonical ensemble MD 
simulation of small systems. The fluid model under consideration is the truncated and 
shifted Lennard- Jones (t. s. LJ) potential with a cutoff radius of 2.5c [30]. Note that the 
chemical potential supersaturation, i.e., S = exp (j8 [fl — fl a (T)]), deviates considerably 
from the pressure supersaturation p/p a and the density supersaturation p/pa, with 
respect to the saturated vapour pressure p<j(T) and density p"{T) of the bulk, cf. Fig. 1. 
For the saturated chemical potential of the t. s. LJ fluid, a correlation based on previously 
published data [8] gives 

fcW-tf) 1-7106E 1.1514c 2 

In Fig. 2, the chemical potential supersaturation is shown as a function of the vapour 
density determined by GCMD simulation with McDonald's daemon. These values agree 
well with the metastable vapour pressure isotherm of the t. s. LJ fluid obtained by 
canonical ensemble simulation. 



INTERVENTION RATE AND NUCLEATION RATE 

The size evolution of any given cluster can be considered as a random walk over 
the order parameter v, changing only by relatively small amounts Av, usually by the 
absorption or emission of monomers. As discussed by Smoluchowski [31, 32] during 
his scientifically most productive period in L'viv and Krakow, the probabilities for the 
growth and decay transitions are proportional to the respective values of the partition 
function W, resulting in 

and 

, 1 (dW/dv)Av 2 . 
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The probability ^ F (v) that a certain size is eventually reached (at any time during the 
random walk process), given that the current size is v, has the property 

^ F (v) = & )+ {v)& >F (v + Av) + & > -(v)& >F (v-Av). (12) 

By substituting 

&> F (v ± Av) = & F (v) ± ^— Av + ^rf^Av 2 + & (Av 3 ) , (13) 

dv 2dv 2 v ' 

it follows for small Av neglecting terms of third order and beyond, that 
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Using the partition function for the grand canonical ensemble, the derivative of the 
probability is given by 

d^ F 

-^ r = rexp(2/3A^2 v ), (15) 

where F is an integration constant. Obtaining the two remaining parameters from the 
boundary conditions 

q\ = 0, (16) 
lim q® = 1, (17) 
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FIGURE 2. Density dependence of the chemical potential supersaturation for the vapour of the t. s. LJ 
fluid, obtained from GCMD simulation with McDonald's daemon (□) and by integration of the Gibbs- 
Duhem equation using data from canonical ensemble MD simulation with T = 0.7 ( — ) and 0.85 e/k^ 
(-)• 



the probability q® for a cluster containing & molecules of eventually reaching macro- 
scopic size, i.e., J? — > °°, is 



/®exp(2j3Afl v )</v 



(18) 
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The intervention rate ^® of McDonald's daemon is related to the nucleation rate Jf by 
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Thus, with an intervention threshold far below the critical size, the intervention rate is 
many orders of magnitude higher the steady-state nucleation rate. However, as confirmed 
by the present simulation results shown in Tab. 1, it reaches a plateau for > v*, where 
v* = 41 according to CNT and 39 according to SPC. 



RESULTS AND DISCUSSION 

Homogeneous nucleation of the t. s. LJ fluid was studied by a series of GCMD simula- 
tions with McDonald's daemon for systems containing up to 17 million particles. 

After a temporal delay, depending on the threshold size, the pressure and the inter- 
vention rate reached a constant value, cf. Fig. 3. In a canonical ensemble MD simula- 
tion under similar conditions as the GCMD simulation that is also shown in Fig. 3, the 
pressure supersaturation decreased from about 3 to 1.5 and the rate of formation was 
significantly lower for larger nuclei, due to the free energy effect accounted for by Eqs. 
(18) and (19) as well as the depletion of the vapour [19]. 
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FIGURE 3. Top: Number per unit volume p n of clusters containing more than 25 (• - ■), 50 ( — ), and 
150 ( — ) particles in a canonical ensemble MD simulation at T = 0.7 e/k$ and p = 0.004044 c~ 3 using 
a hybrid geometric-energetic cluster criterion, number per unit volume p n of clusters with V > 25 (□) in 
a GCMD simulation with T = 0.7 e/k B , S = 2.8658, and = 50, using the Stillinger [22] cluster criterion 
with clusters determined as biconnected components, as well as the aggregated number of McDonald's 
daemon interventions per unit volume in the GCMD simulation, over simulation time. Bottom: Pressure 
over simulation time for the canonical ensemble MD simulation ( — ) and the GCMD simulation with 
McDonald's daemon ( — ) [19]. 



The constant supers aturation of the GCMD simulation agreed approximately with the 
time-dependent supersaturation in the canonical ensemble about t = 400 after simulation 
onset, cf. Fig. 3. At this stage, the number of small clusters present per volume was 
similar in both cases, and the rate of formation for clusters with v > 150 at t = 400 in the 
canonical ensemble simulation was of the same order of magnitude as the intervention 
rate of the daemon. 

TABLE 1. Dependence of the intervention rate as well as the probability 
q@ according to CNT and LFK on the intervention threshold size for McDon- 
ald's daemon during GCMD simulation at T = 0.7 e/kg and S = 2.4958, where 
the rates are given in units of (£/m) 1/2 C7~ 4 . The number of particles in the sys- 
tem and the values for the pressure supersaturation p/p a refer to the steady state 
and the constant volume of the system is given in units of a 3 . 
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Van Meel et al. [16] determined by MC simulation with forward flux sampling that 
supersaturated vapours of the t. s. LJ fluid at a temperature of T = 0.45 e/k^, i.e., 
significantly below the triple point T3 = 0.65 e/k^, initially undergo vapour to liquid 
nucleation, and CNT is known to underestimate the vapour to liquid nucleation rate of 
unpolar fluids [13]. The present daemon intervention rates confirm this conclusion. LFK 
and HSL are significantly more accurate than CNT. Note that in Tab. 2, the nucleation 
rate according to Eq. (19) based on the CNT value of q@ is given. 

From Tab. 2 it is also confirmed that the 'direct observation method' (DOM) [17], 
which in the present case corresponds to assuming 

In J Q = In J - \xiq & = - In xV, (20) 

where x is the temporal delay of formation for the first sufficiently large cluster, is 
inadequate for nucleation near the spinodal line. 

CONCLUSION 

GCMD with McDonald's daemon was established as a method for steady-state simu- 
lation of nucleating vapours at high supersaturations. A series of simulations was con- 
ducted for the t. s. LJ fluid. CNT was found to underpredict the nucleation rate below the 
triple point, whereas LFK and HSL more accurately describe vapour to liquid nucleation 
of the t. s. LJ fluid. 
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TABLE 2. Vapour to liquid nucleation rate at T = 0.45 e/ks from GCMD simulation with 
McDonald's daemon. The theories were evaluated with respect to the metastable vapour-liquid 
equilibrium at p a = 4.28 x 10~ 5 e/a 3 [16], and the vapour-liquid surface tension 7= 1.07 e/a 2 
[16] was used. 
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